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Рассматривается численное решение прямым методом вариационной задачи о построении сверхзвукового сопла с равномер- 
ным потоком на выходе. Предложен способ выбора минимизируемого функционала. Проведено сравнение с результатом, по- 
лученным другим методом. 


1. Введение 

Известно, что если возможна численная (про- 
граммная) реализация некоторой задачи, то эта ре- 
ализация может быть сделана неединственным 
способом. В этом смысле не являются исключени- 
ем и прямые численные методы решения вариаци- 
онных задач, когда решение получается в результа- 
те минимизации определенного функционала. При 
этом для класса вариационных задач, к числу кото- 
рых относятся и газодинамические задачи, некото- 
рые численные реализации могут приводить к то- 
му, что область определения будет являться нео- 
дносвязным множеством. 

В качестве характерного примера рассмотрим 
задачу, решаемую прямым методом, о численном 
построении сверхзвукового сопла с равномерным 
потоком на выходе. Отметим, что к этому близка 
задача о сопле максимальной тяги [1]. Пусть в не- 
котором сопле, профиль которого определяется 
варьируемыми переменными, рассчитывается поле 
течения, а по найденному полю вычисляется функ- 
ционал, характеризующий неравномерность пото- 
ка на выходе из сопла и имеющий минимальное 
значение тогда, когда поток равномерный. 

В основу численной реализации указанной зада- 
чи могут быть положены два основных подхода. 
Первый использует тот факт, что поток остается 


сверхзвуковым и возможно использование простых 
в реализации и быстрых в расчете маршевых схем; 
второй подход допускает существование дозвуко- 
вых течений и, соответственно, требуется примене- 
ние гораздо более сложных и более медленных чи- 
сленных методов. Пусть в основу первой численной 
реализации положен метод характеристик [2]. Тогда 
очевидно, что если в некоторых контурах не будет 
реализовываться полностью сверхзвуковое течение, 
то это приведет к аварийному останову (авосту). 
Следствием этого будет неодносвязность области 
определения, поскольку функционал, использую- 
щий параметры потока на выходе из сопла, вычи- 
слить нельзя. Второй подход, учитывающий по- 
явление дозвукового течения и использующий, на- 
пример, метод установления и схему Годунова [3], 
лишен это недостатка, однако интеллектуальные и 
временные затраты по сравнению с методом харак- 
теристик возрастают на порядок. 

В работе [4] наряду с методом характеристик, 
позволяющим быстро получать решение, предло- 
жен функционал, который, в частности, использу- 
ет значения параметров течения, найденные на 
каждой характеристике С. Это позволяет даже при 
наличии авостов продолжить функционал на одно- 
связную область. 

В данной статье, которая является развитием 
работы [4], предлагается вычислять функционал 
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более простым способом. Для оценки правильно- 
сти полученных результатов использовалось реше- 
ние, полученное по методике [2]. 

2. Постановка задачи 


Необходимо найти контур сопла (рис. 1), на вы- 
ходе из которого в сечении х=х ь поток должен 
иметь однородные поперек сопла параметры 0=0 и 
число Маха М 0 > I . Считается, что давление на вы- 
ходе из сопла больше давления окружающей среды. 


Дано установившееся незакрученное изоэнтро- 
пическое и изоэнергетическое течение идеального 
совершенного газа в осесимметричном сопле, 
рис. 1. Характеристические уравнения и условия 
совместности имеют вид: 


г/х 


= 1 ё(Ѳ±а), 


сІѲ ± 


со^(о) <(о±! іп(а)5ш(ву л=0 . 


(у +1)/2 — со5 2 (а) усо5(0±а) 


где у - показатель адиабаты, далее у=1,4; х и у - 
продольная и поперечная координаты, отнесенные 
к радиусу сопла в минимальном сечении; 
а=аг8Іп( 1/М) - угол Маха, М - число Маха; 
0 - угол наклона вектора скорости (линии тока) к 
оси х; знак +(-) соответствует характеристике 
СДС - ); полагается, что во входном сечении х„=0. 
На стенке реализуется условие непротекания: 
1§(0)=/'(х), где/(х) - функция, описывающая кон- 
тур сопла, а штрих обозначает производную по х. 
Во входном сечении х а поток однородный: 0„=О, 
М т =\ , где индекс «іп» соответствует входу в сопло. 
Далее индекс «0» соответствует выходу из сопла. 



Рис. 1. Схема сопла: оа - минимальное сечение; аЬ - сопло 
с угловой точкой; асі - характеристика С~, принадле- 
жащая пучку волн разрежения; ае - начальная харак- 
теристика ; кі - характеристика СІ; пт - характери- 
стика См; сЬ ~ конечная характеристика С 


При вычислении параметров течения в угловой 
точке «а» используется формула Прандтля-Май- 
ера, следующая из (1): 

Ѳ 2 + у/(а 2 ) = Ѳ } 


у/(а) = -а+ А — — - агсіе 

Л у — 1 


у + 1 

у - 1 


*ё(а) 


( 2 ) 


где нижние индексы 1 и 2 соответствуют параметрам 
до и после разворота в угловой точке; 0,=О, а { =а ы . 


3. Тестовая задача 

При заданном показателе адиабаты рассмотрим 
однопараметрическое, зависящее от М ь семейство 
сопел с угловой точкой и равномерной характери- 
стикой на выходе. Известно [1], что эти сопла име- 
ют максимальную тягу, являются кратчайшими и 
обеспечивают нулевые потери на рассеяние в вы- 
ходном сечении. Поскольку на выходе поток будет 
однородный и параллельный, то указанные сопла 
будут удовлетворять поставленным выше условиям. 

Использование методики [2] при заданных М 0 > 1 
и показателе адиабаты позволяет найти единствен- 
ное сопло, принадлежащее указанному семейству и 
удовлетворяющее условиям задачи. Контур сопла, 
полученный таким образом, считается эталоном, и 
с ним будет проводиться сравнение прямого мето- 
да, а координаты начальной точки «а» и найденной 
конечной точки « Ъ » будут являться исходными дан- 
ными для прямого метода. 

4. Прямой метод 

Даны координаты точек «а», « Ь » и тангенс угла 
наклона сопла в точке «/»>: 

/(Х а ) = Уа » Я х ь) = Уь . І\х ь ) = у'ь, 

*.=0, У а = 1, УІ= 0. (3) 

Последнее равенство следует из условия 0=0 на 
выходе из сопла. 

Для аппроксимации искомого контура в каче- 
стве базисных функций используются степенные 
полиномы: 

т=±сГ, / = (2 х-х ь -х а )/(хь -х а ), 

І=1 

I е [—1,1], хе[х а ,х 6 \. (4) 

Так как должны выполняться условия (3), то ко- 
эффициенты с т , Сдг+ 2 , с т и профиль сопла Лх) вы- 
ражаются через N линейно независимых коэффи- 
циентов с ь ..., с, ѵ . Требуется найти профиль сопла, 
удовлетворяющий геометрическим условиям (3) и 
постановке задачи. 

Прямой расчет. В рамках прямого метода приня- 
то называть прямым расчетом единичный расчет 
поля течения, который заканчивается вычислени- 
ем функционала. Выбор функционала будет обсуж- 
даться ниже. В прямом расчете система уравнений 
(1) решается по схеме, рис. 2, где цифрами «1» и «2» 
обозначены точки с известными параметрами, из 
которых выходят характеристики С* и С . Цифрой 
«3» обозначена точка их пересечения. Записывая 
уравнения ( 1) в разностном виде, получим систему 
уравнений относительно неизвестных параметров 
а/, 0/, хі, у{ в точке «3»: 
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Уз -У 2 


Уз -У, 


= іё(Ѳ 23 -а 23 ), 


хі -х. 


= 1§(0 13 +а 13 ), 


0 /- 0 ,+ 


со з а, 


(у + 1)/2-со8 2 а 
8Іпа,, зіп 0, 


-(а/ -а,) + 


Уіз со8(Ѳ 13 +а 13 ) 


(х 3 ~х,) = О, 


ѳі-ѳ 2 - 


со8~ а. 


(7 + 1)/2-со 8 2 а 2 


-(а/ -а 2 )- 


8Іпа„8ІпѲ^ . : . л 

- (х 3 —х 2 ) = О, 


[ 3^ 23 С08 (^23 ^ 23 ) 


здесь у=1,2,... - номер итерации. Обозначим 
р=а,Ѳ,х,у, тогда р 13 =(р | +р/-‘)/2, р а =(р 2 +рПЯ По- 
лученная система решалась итерационно до вы- 
полнения условия тах|р/-/?/ _1 |<1(Ь 8 . На начальной 
итерации полагалось Л 0= (Рі+Л)/2- Так как на оси 
0=0, у=0, а на стенке 0=агс1§(/' , (л;)) и у=/{х), то в 
данной системе уравнений производились очевид- 
ные упрощения. 



Рис. 2. Схема расчета: 13 ~ отрезок характеристики С; 
23 - отрезок характеристики С 


Таким образом, от оси к стенке, по известной ха- 
рактеристике С* рассчитывается следующая харак- 
теристика С? +1 , где индексы і и /+1 — порядковые но- 
мера вычисляемых характеристик. Расчет произво- 
дится до тех пор, пока очередная характеристика С 
не приходит в точку «Ь». В качестве начальных усло- 
вий используются 0 Й =О и М іп = 1,001 на начальной ха- 
рактеристике ае. Такие начальные условия достаточ- 
но часто применяют вместо плоской звуковой линии 
[5]. Отметим, что для вычисления интегралов ис- 
пользовалась формула трапеции. 

Выбор функционала. В [4] использовался функ- 
ционал, который рассчитывался вдоль характери- 
стики (рис. 1): 



( 5 ) 


Если при вычислении характеристики С* +1 проис- 
ходил авост, то формула (5) заменялась выражением: 


■! =(х„-х,)+ 



<р(х) = 


Г|/(*)|, /'(*)< о 

|о, /'(*)> 0 


( 6 ) 


где первый интеграл вычислялся вдоль последней 
рассчитанной характеристики С, + , а второй - сум- 
мировал площадь участков контура сопла с отрица- 
тельными углами наклона. Если авост произошел 
при расчете характеристики С[ + , то х,=0 и 0=0 вдоль 
начальной характеристики С 2 . В результате функ- 
ционал (6) приобретает вид: 

ь 

/ = х ь + ^(р(х)сіх. (7) 

а 

Как видно, способ вычисления функционала 
(5-7), предложенный в [4], достаточно сложен. И, 
кроме этого, замена формулы (5) на формулу (6) 
приводит к тому, что функционал, сохраняя непре- 
рывность, в общем случае может стать недиффе- 
ренцируемым. 

В данной работе предлагается другой подход. 
Рассмотрим два функционала 



которые вычисляются вдоль произвольной харак- 
теристики С*. Выбор функционалов (8) обусловлен 
следующими соображениями. 

Пусть сопло с известным М 0 удовлетворяет 
условиям задачи. Тогда, как видно на рис. 3, пер- 
вый функционал из (8) имеет два минимума /[=0. 
Первому минимуму на характеристике ае соответ- 
ствует 0=0 и о: ( „=агс8т(1/44). Второму минимуму 
на характеристике сЬ соответствует 0=0 и 
а 0 =аг8Іп(1Д/ 0 ). Как следует из уравнений (1), толь- 
ко второй минимум является необходимым и до- 
статочным условием решения задачи. Очевидно, 
что использование первого функционала не гаран- 
тирует единственности решения. 



Рис. 3. Значение функционалов вдоль характеристик С*, 
і - номер характеристики 

Второй функционал из (8), в отличие от перво- 
го, имеет единственный минимум /,=0 на характе- 
ристике сЬ (рис. 3), чему соответствует а=а 0 . Одна- 
ко этот минимум является необходимым, но не до- 
статочным условием решения задачи, поскольку 
может не выполняться условие 0=0. Комбинация /, 
и / 2 дает следующий функционал: 
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который имеет единственный минимум / 3 =0 на ха- 
рактеристике сЪ (рис. 3). Так как этот минимум яв- 
ляется необходимым и достаточным условием ре- 
шения задачи и ему соответствует Ѳ=0 и а=а 0 , то 
далее будет использоваться функционал (9). 

При вычислении функционала (9) воспользу- 
емся приемом, предложенным в [4], и если при вы- 
числении характеристики С* +1 происходит авост, то 
функционал (9) вычисляется вдоль уже рассчитан- 
ной характеристики Ц. 

Если сопло имеет отрицательный наклон сопла 
/\х а )<0 в точке «а», то авост происходит при вычи- 
слении характеристики С ъ и функционал (9) вычи- 
сляется на начальной характеристике ае. Посколь- 
ку при вычислении параметров течения в угловой 
точке «а» используется формула Прандтля- Майера 
(2), то можно использовать эту формулу и для вы- 
числения а 2 при Ѳ 2 <0. В этом случае найденные 
значения а 2 не имеют физического смысла, т. к. 
этому соответствует дозвуковое течение. Однако в 
силу своей непрерывности (рис. 4) функция а 2 (Ѳ 2 ) 
позволяет продолжить функционал (9) на односвя- 
зную область в случае отрицательного наклона со- 
пла /'(х о )<0 и обеспечить функционалу непрерыв- 
ность и дифференцируемость. Как показали расче- 
ты, при вычислении по формуле трапеции функ- 
ционала на начальной характеристике ае достаточ- 
но брать только две точки - «а» и «е». 



Поскольку на входе в сопло граничные условия 
не меняются, а параметры течения определяются 
только профилем сопла (4), то значение функциона- 
ла (9) будет неявно зависеть от этого профиля. Таким 
образом, задача нахождения профиля сопла /(х), до- 
ставляющего экстремум функционалу (9), сводится 
К поиску ТОЧКИ (с ь ...,Сц), в которой функция многих 
переменных /=/,((/,. ..,с іѴ ) имеет экстремум. Для на- 
хождения минимума этой функции использовался 
квазиньютоновский метод Бройдена из [6]. 

5. Численные результаты 

Число N варьируемых переменных с, менялось 
от 1 до 10, начальное значение с,=0. На характери- 


стике ае задавалось 50 точек. Координаты точек «а» 
и «Ъ»\ (0; 1) и (3,576; 1,299). Этим значениям соот- 
ветствует найденное из [2] эталонное сопло с угло- 
вой точкой и числом Маха на выходе М (] =2. 

Сравнение решения для обоих функционалов с 
эталонным соплом показало, что при 7Ѵ= 1 0 орди- 
наты сопел отличались в четвертом знаке после за- 
пятой, а максимальная относительная погреш- 
ность по составила 0,02 %. 

На рис. 5 представлены результаты минимиза- 
ции в виде окончательных значений функционалов 
для различного числа N. Таблица демонстрирует 
влияние числа коэффициентов и используемого 
функционала как на число прямых расчетов, так и 
на количество авостов. Из представленных резуль- 
татов видно, что предложенный функционал (9), 
несмотря на простоту, по эффективности не усту- 
пает функционалу (5-7), предложенному в [4]. 
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Рис. 5. Минимизация функционалов: 1) формулы (5-7) из 
[4], 2) предлагаемая формула (9) 


Таблица. Минимизация функционалов. В числителе ~ чи- 
сло прямых расчетов, в знаменателе - авостов 


№ 

Формулы (5-7) [4] 

Предлагаемая формула (9) 

1 

13/2 

17/2 

2 

36/3 

39/4 

3 

63/4 

69/5 

4 

98/4 

99/5 

Б 

135/Б 

115/7 

6 

150/7 

166/7 

7 

199/7 

203/6 

8 

273/6 

304/7 

9 

385/8 

320/8 

10 

444/7 

492/9 


Заключение 

Проведенные численные исследования показы- 
вают, что предложенный функционал достаточно 
прост в реализации, при наличии авостов имеет 
продолжение на односвязную область и сохраняет 
дифференцируемость, а по эффективности не усту- 
пает ранее рассмотренному в работе [4] функцио- 
налу. Сравнение профиля найденного сопла с эта- 
лонным профилем показало, что при 10 варьируе- 
мых переменных максимальная относительная по- 
грешность по ординате составляет 0,02 %. 
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